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Abstract — The capacity or approximations to capacity of var- 
ious single-source single-destination relay network models has 
been characterized in terms of the cut-set upper bound. In 
principle, a direct computation of this bound requires evaluating 
the cut capacity over exponentially many cuts. We show that the 
minimum cut capacity of a relay network under some special 
assumptions can be cast as a minimization of a submodular 
function, and as a result, can be computed efficiently. We use 
this result to show that the capacity, or an approximation to the 
capacity within a constant gap for the Gaussian, wireless erasure, 
and Avestimehr-Diggavi-Tse deterministic relay network models 
can be computed in polynomial time. We present some empirical 
results showing that computing constant-gap approximations to 
the capacity of Gaussian relay networks with around 300 nodes 
can be done in order of minutes. 

For Gaussian networks, cut-set capacities are also functions 
of the powers assigned to the nodes. We consider a family of 
power optimization problems and show that they can be solved 
in polynomial time. In particular, we show that the minimization 
of the sum of powers assigned to the nodes subject to a minimum 
rate constraint (measured in terms of cut-set bounds) can be 
computed in polynomial time. We propose an heuristic algorithm 
to solve this problem and measure its performance through 
simulations on random Gaussian networks. We observe that in 
the optimal allocations most of the power is assigned to a small 
subset of relays, which suggests that network simplification may 
be possible without excessive performance degradation. 

Index Terms — capacity, network simplification, power alloca- 
tion, relay networks, submodular optimization. 



I. Introduction 

Relay networks, where one or more source nodes send 
information to one or more destination nodes with the help 
of intermediate nodes acting as relays, are often used to 
model communication in wireless sensor networks. In sensor 
networks, sensor nodes have limited power sources and often 
require multi-hop communication with the help of intermediate 
nodes to reach the data aggregation centers. To guide the 
design of these networks it is of interest to characterize 
fundamental communication limits such as the capacity, which 
represents the maximum reliable communication rate. 

Various communication models for relay networks capture 
in an abstract setting different aspects of practical systems. 
The wireless erasure network model of [8 1 captures the effect 
of packet losses in the wireless setting. The deterministic 
network model of Avestimehr, Diggavi and Tse (ADT) J4) 
incorporates broadcast and interference and can be used to 
gain insights about communication in more complex models 
that incorporate noise. Among these, of special importance 
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is the Gaussian relay network, which models power limited 
transmitters and received signals corrupted by additive white 
Gaussian noise. 

While the capacity of some network models (e.g. wireless 
erasure and ADT) is well characterized, the capacity of the 
Gaussian relay network, even in its simplest form with one 
transmitter, one relay, and one receiver, is in general unknown. 
The best known capacity upper bound is the so-called cut-set 
bound. A cut O of a network can be considered as a subset 
of nodes which includes the source node and excludes the 
destination node. For this cut, the capacity F(il) is defined 
as the maximum rate that information can be transferred form 
the nodes in f2 to the nodes that are not in conditioned on 
the fact the information on il c (the nodes that are not in ft) is 
known. The cut-set upper bound is the minimum cut capacity 
over all the possible cuts. 

In the Gaussian setting, there are several capacity lower 
bounds based on different communication schemes, such 
as amplify-and-forward, decode-and-forward, compress-and- 
forward, quantize-and-forward, etc. @, Q, lfl8l . Recently, 
Avestimehr, et al. ||2) made significant progress in the capacity 
characterization of Gaussian relay networks by showing that a 
quantization and coding communication scheme can achieve a 
communication rate within a constant gap of the cut-set upper 
bound, where the gap only depends on the number of nodes 
in the network (i.e. it is independent of the channel gains 
and power levels). However, the evaluation of the achievable 
communication rate, which is necessary to implement the 
scheme, requires the computation of the cut-set bound for the 
network. Assuming that for a given cut the cut capacity is 
easy to compute, finding the cut-set upper bound can be a 
challenging problem. For a network with n relays there are 
2™ different cuts and a greedy algorithm needs exponential 
time in the number of relays to compute the cut-set bound. 

In this work we show that the achievable rate of the scheme 
of Q for the Gaussian relay network can be computed in 
polynomial time, and as a result, can be computed efficiently. 
This result is obtained by showing that the cut capacity of a 
fairly large class of networks under the assumption of inde- 
pendent encoding at the nodes in fl is a submodular function. 
For the special case of layered relay networks, ll27ll showed 
the equivalent of our submodularity result simultaneously with 
our conference version of this paper l25l . Submodularity 
properties of conditional entropy (in terms of which cut- 
capacities are expressed) have also been used in 1211, 1131 to 
bound the cut-capacity of a network in terms of the cut- 
capacity of the corresponding unfolded graplfl 

'Please, refer to (2| for the definition of an unfolded graph. 
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Existing results on minimization of submodular functions 
provide algorithms with polynomial running time 0(n 5 a + 
n 6 ), where a is the time that it takes to compute F(Q) and 
n is the number of nodes in the network ll20l . In addition, 
there exist possibly faster algorithms without polynomial time 
performance guarantees based on Wolfe's minimization norm 
algorithm [12|. In Section [VTl by simulations, we show that 
the cut-set bound for a Gaussian relay network with around 
300 nodes can be computed on a laptop computer in about a 
minute using a Matlab package for submodular minimization 
provided in |[T9l . 

Our results, extend and generalize previous results for the 
ADT model. This model can be seen as a high signal- 
to-noise-ratio (SNR) approximation of the Gaussian model, 
incorporating the effects of broadcasting and superposition of 
signals while de-emphasizing the effects of noise. Amaudruz 
et al. HI showed that the cut-set bound for a layered ADT 
model can be computed efficiently. They have extended graph 
flow algorithms such as Ford-Fulkerson's in a nontrivial way 
to find the maximum possible linearly independent (LI) paths 
in the network. They showed that the capacity of the network 
is equal to the maximum number of (LI) paths and can 
be computed in time 0(M ■ \E\ ■ C 5 ), where M is the 
maximum number of nodes per layer, \E\ is the total number 
of edges and C is the capacity of the network. Moreover, 
they showed that the capacity can be achieved by using 
a very simple one-bit processing at the relay nodes. Later 
Goemans et al. |[T3l showed that the deterministic model is 
a special case of a flow model based on linking systems, a 
combinatorial structure with a tight connection to matroids. 
As a by-product, they obtained the submodularity of the cut 
capacity for layered ADT networks. Using this observation 
they provided various algorithms related to matroid theory to 
compute the cut capacity of the layered deterministic model 
based on finding intersection or partition of matroids. These 
results led to faster algorithms to compute the capacity of 
large layered ADT networks. In addition, there has been 
other extensions on improving the running time of the current 
algorithms for computing the capacity of ADT networks (9), 

ed, ma, ma. 

In addition to showing that the capacity within a constant 
gap of the Gaussian relay network can be computed in polyno- 
mial time, our results allow us to compute in polynomial time 
the capacity of the wireless erasure network. Furthermore, we 
provide a simple proof for the computability in polynomial 
time of the capacity of the layered and non-layered ADT 
networks. 

Building on the submodularity of the cut-capacity for 
independent encoding at the nodes, we show that, in the 
Gaussian setting, it is possible to efficiently optimize the power 
allocated to the source and relay nodes. We consider two 
power optimization problems: (i) minimize the total power 
satisfying a minimum source-destination data rate constraint 
and power constraints at each node; (ii) maximize the source- 
destination data rate satisfying total and individual power 

2 In a layered network, the nodes in one layer are only connected to the 
nodes in the next adjacent layer. In particular, there is no direct connection 
from source to destination. 



constraints at the nodes. Since the capacity of the Gaussian 
relay network is approximately given by the cut-set upper 
bound with independent encoding at the nodes, we use this 
cut-set bound to characterize data rate in the optimization 
problems. We show that these optimization problems can be 
solved in polynomial time and use simulations to get insights 
about some properties of the optimal power allocations for 
networks of various sizes. We observe that optimal power 
allocations assign most of the power to a small subset of 
nodes and that setting the power to zero in the remaining nodes 
(i.e. removing these nodes from the network) often results in 
a small rate loss. Nazaroglu, et al. showed in [26 1 that for 
the special case of the iV-relay Gaussian diamond network a 
fraction k/(k + 1) of the total capacity can be approximately 
achieved by using only k of the total N relays. This suggests 
that the diamond network can be significantly simplified by 
tolerating a small performance loss. Our results provide a 
numerical counterpart to the fundamental performance bounds 
derived in [26 1 and suggest that network simplification may 
also be possible in more general Gaussian relay networks. 

We obtain these results by considering a general framework 
to compute the cut-set bound. We assign transmit signal ran- 
dom variable Xi to node i £ V and we assume the probability 
distribution over the signals Xi,X 2 , ■ ■ ■ ,X n to be indepen- 
dent, i.e p(X 1 ,X 2 , ...,X n )= p 1 {X 1 )p 2 {X 2 ) ■ ■■p n (X n ). We 
also assign received signal random variables Yj's to each 
node. The network is defined by the transition probability 
function /(Yi, Y 2 , Y n \X 1 ,X 2 , ■ ■ ■ , X n ). We further as- 
sume that the transition probability function is of the form 
fi{Yx\X u X n ) ■ ■ ■ f n (Y n \X u X n ), meaning that the 
received signals are independent conditioned on the trans- 
mitted signals in the network. For such networks we show 
that F(Q) = /(Yfjc; Xn|Xocj^| is submodular with respect 
to fi. Later we show that for ADT networks, the Gaussian 
relay network and the wireless erasure network, we can find 
Px(Xi) ■ ■ ■ p n (X n ) such that minn F(Q) becomes equal to 
the capacity or the capacity within a constant gap. In other 
words, the min-cut problem for these networks can be cast as 
a minimization of a submodular function. 

The paper is organized as follows. In Section Hill we show 
that for specific type of networks the cut value, F(ft), is a 
submodular function. We then show in Section [TV] that for 
many wireless network models such as the ADT deterministic 
network, Gaussian relay network and wireless erasure network 
the capacity or an approximation to the capacity can be 
cast as a minimization of F(H,). In Section [V] we study 
two power optimization problems and show that they can be 
solved efficiently. Finally, in Section |VI] we describe results 
related to solving optimization problems involving submodular 
functions and perform power optimization in various randomly 
generated networks of different sizes. We start by introducing 
the notation used in the rest of the paper. 

II. Notation 

Let V denote the set of nodes in the network and |V| its 
cardinality. For any subset A of nodes we denote by V\A 

3 See Section [ii] for a definition of the notation Xq, Yqc, etc. 
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or A c the set of nodes in V that are not in A. We assume 
V\A U B = V\(A U B). A cut O is defined as a subset of 
nodes in V. A cut splits the nodes in the network into two 
groups, the nodes that are in fl and the ones that belong to 
V\f2. Random variables are shown in capital letters such as Xi 
and Yi. We use boldface letter for vectors, e.g. x is a constant 
vector and X is a random vector. We use Xji to denote 
(X Vl ,X V2 ,..., X„ m ) with Vi G ft. The function I(X; Y\Z) 
is the mutual information between random variables X and 
Y conditioned on random variable Z. With a slight abuse 
of notation we use H(X) to denote either the entropy or 
differential entropy of the discrete or continuous random 
var iableX @. By ¥ p we denote a finite field with p elements. 
Finally, all the log(-) functions are in base two. 

III. SUBMODULARITY OF CUT-SET FUNCTION 

Submodularity arises in many combinatorial optimization 
problems and large body of research has been developed 
on minimizing or maximizing submodular functions under 
various constraints. 

A submodular function / : 2 V — > R is defined as a function 
over subsets of V with diminishing marginal returns, i.e. if 
A,B CV with A C B and any v G V\B, 

f(AUv)-f(A)>f(BUv)-f(B). 

The theorem below establishes the submodularity of the 
cut capacity function of a general relay network under some 
special assumptions. This theorem will be used in Section [IV] 
to prove that the capacity or an approximation to the capacity 
of various specific relay network models can be computed by 
minimizing a submodular function. 

Theorem 1. Consider a network consisting of nodes in V. 
Each node sends a message Xi,i G V and receives Yi,i G 
V. If the messages are independent p{X\, X2, X\ v \) = 
Pi(Xi)p2(X2) • ■ • p\v\ (^|V|) and conditioned on the sent mes- 
sages the received messages are independent, then the function 

F(A) = I(X A ;Y V \ A \X V \ A ) , AC V 

is submodular. 

Proof: To show that F(A) is submodular we show that 
F(A U a) — F(A) is monotonically non-increasing in A for 

a <£ A. 

F(AUa) = I(X Aua ;Y v \ AlJa \X v \ Aua ) 

^= H(X. Aua \~K v \ Aua ) — H (X. Aua \Y v \ Alla , X v \ Aua ) 

= HQLa) + H(X a \~X.A) — H(X a \Y V \ Aua , X V \AUa) 

- H(X A \X a , Yyy^ijQ, Xyy^ua) 

=H(X A ) + H(X a \X A ) - i^XalYvytu^Xv^ua) 

- H(X A \Y v \ Aua , Xyyt) 

where (a) is the definition of mutual information and (b) is 
from the chain rule for the entropy function. Therefore, 

F(Al)a) - F(A) 

=H(X a \'X A ) - H(X a \Y v \ Aua , ~K-v\Aua) 



- H(X A \Y v \ Aua , X. V \ A ) 

+ H(X A \Y v \ AUa ,Y a ,X v \ A ) 
=H(X a \X. A ) - H(X a \Y v \ AUa , X v \Aua) 

- I(X. A ; Y a \Y v \ AUa , X V \ A ) 
=H(X a \X. A ) - H(X a \Y v \ Aua , X v ^ Ua ) 

- H(Y a \Y v ^ AUa , X V \ A ) 

+ H (Y a \X A ,Y v \ Alla ,X v \ A ) 
= J?(X„|X^) — H(X a \Y v \ AUa , Xv\a) 
non-increasing in a nondecreasing in A 

-H(Y a \Y v \AUa,X V \A)+H(Y a \X V ) 

nondecreasing in A 

where the last equality follows because Y a is independent of 
Yy\Aua conditioned on Xy. So, F(A U a) — F(A) is non- 
increasing in A and thus F(A) is submodular. ■ 
In the following example we show that if the signals 
at the nodes are correlated then F(A) is not necessarily a 
submodular function. 

Example. Consider a symmetric Gaussian diamond network 
with two relays such that the channel gains from source to 
relays are equal to one and from relays to destination are equal 
to three. Letting X s ,X ri , and X r2 be the signals transmitted 
at the source and relay nodes, then the received signals at 
relays and destination are given by 

Yr, \ / 1 \ / X, \ / Z„ \ 

Y- r i 3 3 X ' Hz'* 
id/ \ U 6 6 J \ A r2 J \ Zi d J 

where Z ri ,Z r2 ,Z ( j are i.i.d. A/"(0, 1). For this example, we 
set the probability distribution of X^X^X,,,, to be jointly 
Gaussian with zero mean and covariance matrix 

/I 0\ 
S= 1 p ]. 

\ P 1 / 

Finally, consider the sets A = {s, r\] and B = {s,r2}. 
Figure [T] shows how the function F(A) + F(B) - F(A U 
B) — F(A n B) varies for different values of the corre- 
lation coefficient (between X ri and X r2 ) p G [0, 1]. We 
see that F(A) + F(B) can be greater than or less than 
F(AUB) +F(Ar\B) depending on the value of p. It follows 
that in general F(-) is not a submodular or a supermodular 
function when the there is correlation among the signals at 
the nodes. 



IV. Wireless network models 

In this section, by applying the result of TheoremQ] we 
show that the capacity or an approximation to the capacity for 
the ADT deterministic network, Gaussian relay network, and 
wireless erasure network can be cast as a minimization of a 
submodular function. 
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Fig. 1. F(A)+F(B)-F(AUB)-F(AnB) as a function of the correlation 
coefficient p. In general, F(-) is neither submodular or supermodular when 
the signals are correlated. 

A. Deterministic model (ADT) 

We start by briefly describing the network model of [40. In 
this model, each link from node i to node j has an associated 
non-negative integer gain ny. Each node i £ V transmits a 
signal X, and receives a signal Y,, both in where q = 
maxjj riij. At any given time, the received signal at node j 
is given by 

Y, Y. S ' ? " X < 

i£V\{d} 

where d is the destination node, the shifting matrix S is given 
by 

/ •■■ \ 
10 ••• 
s= 1 ••• 

\ ■•■ 1 / 

and the sums and products are in F p . 

For a given cut SI of the network, where SI includes the 
source node and excludes the destination node, we can stack 
together the input vectors , i G SI and output vectors 
Yi,i £ Sl c , and define a transition matrix Aq that gives the 
input-output relationship of these vectors according to (Q~|). It 
is shown in [2 | that the capacity of the deterministic network 
is equal to miriQ rank(An). We show next in Theorem [2] 
that rank(Ao) is submodular, and hence the capacity can be 
computed by minimizing a submodular function. 

Proposition 1. Assume an m x n matrix A over ¥ p . Let Af 

be the subspace Af == {x £ F p ' | Ax = 0}, and let be the 
set of cosets of Af in F™. Pick to be an element in the ith 
coset of Af for i = 1, 2, . . . , \&\, and set = Axj. Notice that 
Yi Yj if i 3- Now, if we choose x uniformly at random 
from elements o/F™ with probability 1/|F™|, then the mapping 

4 Please, refer [4] for a more complete description of the model and its 
motivation. 



Ax maps x to {yi,y2, • • ■ uniformly at random with 

probability l/|Sf|. In addition, the cosets of Af form a partition 
o/F£ into p n /\Af \ sets. Also rank(A) +log p (\Af\) = n. Thus, 
log p = rank(A). 

Theorem 2. For a deterministic model, given a cut SI assume 
Aq is the transition matrix form nodes in fl to nodes in f2 c . 
Set D(Cl) = rank(Aii), then D(Q) is submodular. 

Remark 1. A special case of Theorem^ for layered ADT 
networks was proved in earlier works l\13V , H23V . 

Proof: In the network, assume node i sends 6j sym- 
bols Xi t i,Xi t 2, . . . ,Xij,t with Xij G F p . We assume 2i,/s 
drawn i.i.d. with uniform probability distribution over F p , 
i.e. p(xij = q) = l/|F p | for all q G F p . From the 
definition of transition matrix, An, if we assume for the cut 
f2, s = (si, S2, ■ ■ ■ , Sfe)* symbols are being sent from nodes 
in n and r = (n, f2, ■ ■ ■ , r^)* symbols are being received by 
nodes in 57 c then r = Aqs. Then we can write 

7(Xo;Yac|Xrj c ) =H(Y n o\Xn°) - ff(Y n =|X n ,X n c) 

^(YnelXn.) 
=^(A n s|Xoc) 

^logp |S? | =rank(A ) 

where 5f is the set of cosets of Af where Af = {s : Aqs = 0}. 
Equality (a) is because Yqc is a deterministic function of Xq 
and (b) is the result of PropositionQ] and the fact the s has 
uniform probability distribution. 

Notice that for the independent probability distribution 
on the sources the received signals are independent con- 
ditioned on transmitted signals so, based on TheoremQ] 
/(Xn; Yfie|Xn<=) which is equal to D(il) is submodular. ■ 

B. Gaussian relay network 

The Gaussian network model captures the effects of broad- 
casting, superposition and noise of power constrained wireless 
networks. In this model, at any time index (which we omit) 
the received signal at node j G V\{s} is given by 

)) J2 hjXi + Nj (2) 

i£V\{d} 

where X.- L G C is the transmitted signal at node i, subject to 
an average power constraint E{\Xi\ 2 ) ^ 1, hy G C is the 
channel gain from node i to node j, and Nj G CAf(0, 1) is 
additive white circularly symmetric complex Gaussian noise, 
independent for different j. 

It has been show in ||2T1 Theorem 2.1] that using lattice 
codes for transmission and quantization at the relays, all rates 
R between source {s} and destination {d} satisfying 

H<min/(Xn;Ync|X U c)-|V| (3) 

can be achieved, where SI is a source-destination cut of the 
network and Xq = {X u i G SI} are i.i.d. CAf (0,1). In 
addition, the restriction to i.i.d. Gaussian input distributions is 
within |V| bits/s/Hz of the cut-set upper bound [2|. Therefore 
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the rate achieved using lattice codes in the above result is 
within 2 1 V | bits/s/Hz of the capacity of the network. 

The following corollary is an immediate consequence of 
TheoremQ] 

Corollary 1. The function F(Q) = /(Xq; Yqc |Xn=) with the 
elements o/Xo being i.i.d. £/V(0, 1) is submodular. 

Due to CorollaryO] the minimization in (O is the minimiza- 
tion of a submodular function and the resulting optimal value 
is within 2|V| of the capacity of the network|5 

C. Wireless erasure network 

In (8) the authors introduce a special class of wireless 
networks, called wireless erasure networks. In these networks, 
a directed graph Q — (V, £ ) defines the interconnections 
between nodes. To model the broadcast effect of wireless 
networks, the signals on all outgoing arcs of any given node 
are equal to each other. There is no interference among 
multiple arcs arriving at a given node in this model, and the 
signals on the various arcs are erased independently of each 
other. We assume binary transmitted signals at each node, i.e. 
Xi G {0, 1}, i G V\{d}, but all the results can be extended to 
models with larger input alphabets. It has been shown in [8| 
that the capacity of the network is 



G = minF(fi) = min£ 1- J] Hi 



(4) 



where e, j is the probability of erasure when node i is sending 
information to node j. We show in the following theorem that 
F(fl) is submodular. 



Theorem 3. The function F(£l) = J2ien ( 1 - ]ljen<= e *j 
equals /(Xq; Yjje |X^c) where Xi are i.i.d. ~ Bernoulli(\j1) 
for i G ft. Therefore, F(Cl) is submodular. 

Proof: For i.i.d. ~ Bernoulli (1/2), we can write 

^(Xfj;Yoc Xqc) 

=fr(Xn|Xn.)-fr(Xn|Y n .,Xn.) 



®X)(l--ff(^i|Y n .)) 



HWYj = yj,j G n c )p(Yj = yj , j G O c ) 
= 53 (l - HiXilYj = e,j G O c )p(Yj =e,je 



(c) 



Eh n 



5 Notice that Z(X n ; Y n c |X n <=) = logdet(J + ffift) where _ff is the 
matrix of channel gains from nodes in f2 to nodes in Q c and if' is the 
conjugate transpose of H. Therefore, it is easy to compute the capacity of 
each cut. 



We used in (a) the independence among Xi and the channel 
erasures, in (b) the fact that for Xi ~ Bernoulli(l/2), H(Xi) = 
1, and in (c) the fact that for Xi ~ Bernoulli 1/2), H(Xi\Yj = 
e, j G £l c ) = 1 and for independent erasures we have p(Yj = 
e, j G £l c ) = ri^enc e u- Theorem[T]can be applied to conclude 
that F(Q,) is submodular. ■ 



V. Power optimization 

In the previous section, for the Gaussian relay network 
model, we considered fixed power assignments to the different 
nodes in the network, and have shown that a constant gap 
approximation to the capacity can be efficiently computed. 
In many applications it is of interest to allocate the nodes' 
transmission powers to optimize a given objective. For ex- 
ample, in a network where the nodes are battery powered, it 
may be of interest to maximize the network lifetime while 
satisfying a baseline quality of service. Alternatively, it may 
be desirable to maximize the network throughput for a given 
total power budget. This total power budget may arise due to, 
e.g., a maximum system weight constraint which is dominated 
by the battery weight, or a total system cost, which may be 
heavily influenced by the cost of the batteries. Power allocation 
optimization may also naturally arise in situations where the 
channel gains, while known to all the nodes, slowly vary over 
time. In this case, it may be desirable to optimally allocate 
power for the current channel condition. 

As before, we characterize communication rates in terms 
of cut-set capacities. We consider a model where the cut-set 
capacities are functions of the cuts and powers assigned to the 
nodes in the network: F(Q, p) : V x IRl v l — > R, and we focus 
on the Gaussian model where this function depends explicitly 
on the power assignment, 

F(n,p) = J(Xn;Y n .|Xn.) = logdet(/ + fr n flu.H&) (5) 

where Hq is the matrix of channel gains from nodes in £1 to 
nodes in fl c , W is the conjugate transpose of H, and Pq is 
a diagonal matrix where the diagonal elements are the powers 
of the nodes in Q. 

We will show in LemmaQ] below that in the Gaussian case 
F(fl,p) is a concave function of p. While the results of this 
section are stated and proved for the Gaussian model, we 
conjecture that similar results should hold for other models 
in which F(ft, p) is a concave function of p. 

For Gaussian relay networks, we show that the following 
optimization problem can be solved in polynomial time, 



minimize 

R,P,P 

subject to 



R < F(fl U {s}, p) for all Q C V\{d} 

< P < Pmax 
|V| 

i=l 

Ro < R, P ^ P tot 



(6) 



for fixed constants /i 1 ,/^ 2 ,i?o> Pmax an d Ptot- In me rest of 
the section, we denote the feasible set of the optimization (|6]l 
by K. 
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We use the Ellipsoid method |[T4l . ITT31 to show that the 
optimization §6§ can be solved efficiently. We will use the 
following definitions and result. The reader is referred to |fT31 
for more details. 

Definition 1 (Polynomial computability). A family of opti- 
mization programs is polynomially computable if: 

(i) for any instance of the program and any point x in 
the domain, the objective and its subgradient can be 
computed in polynomial time in the size of the instance. 

(ii) for a given measure of infeasibility Infeas(-), it should 
be possible to determine if Infeas{n) e in polynomial 
time, and when this inequality is not satisfied, it should 
be possible to find in polynomial time a vector c such 
that 

c T x > c T y, Vy : Infeas(y) ^ e. 

Definition 2 (Polynomial growth). A family of optimization 
programs has polynomial growth if the objectives and the 
infeasibility measures as functions of points x in the domain 
grow polynomially with ||x||i. 

Definition 3 (Polynomial boundedness of feasible sets). A 
family of optimization programs has polynomially bounded 
feasible sets if the feasible set of an instance of the program 
is contained in an Euclidean ball centered at the origin with 
radius that grows at most polynomially with the size of the 
instance. 

Proposition 2 (JT5] Theorem 5.3.1]). Let V be a family 
of convex optimization programs equipped with infeasibility 
measure Infeas(-). Assume that the family is polynomially 
computable with polynomial growth and with polynomially 
bounded feasible sets. Then V is polynomially solvable. 

In order to use Proposition^ we need to check that the 
optimization © is a convex program. Since the objective 
function is linear, we only need to check that the feasible 
set K. is convex. 

Lemma 1. The feasible set /C is a convex set. 

Proof: First we show that the function P(f2, p) = 
logdet(7 + iJjiPn-fffj) is concave in p where ^ p for any 
cut O C V. For any two vectors pi, p 2 and 7 € [0, 1] we 
can write 

7F(n,pi) + (i-7)f(n,P2) 

= 7 logdet(J + tf n Pii4) 

+ (l- 7 ) logdet(/ + J ffaP 2 -ff+) 

I ]ogdBk(rf(I + HnPtHk) 

= logdet(/ + Hn^Px + (1 - 7)^2)^) 

where Pi (P2) is a diagonal matrix where the diagonal 
elements are the elements of pi (P2) that belong to f2 
(respectively), and (a) follows from the concavity of log det X, 
for X >- ED- 

Next, consider the set = {(P,P,p) : R < F(ft,p)}. 

Choose two vectors (Pi,Pi,pi) and (P 2 ,P 2 ,p 2 ) in ^(O). 



We will show that ^(O) is a convex set by showing that for 
any 7 G [0,1] the vector 7 (Pi, Pi, pi) + (l- 7 )(P 2 ,-P2,P2) is 
also in #(0). Notice that Pi < F(Q, pi) and P 2 < F(fi, p 2 ) 
if and only if (JJi,P ljPl ) G and (P 2 ,P 2 , P2) 6 *?(fi). 

Therefore 

7P1 + (1 - 7 )P 2 < 7 P(fi, Pi) + (1 - j)F{n, p 2 ) 
( | ) P(0,7P1 + (1- 7 ) P2 ) 

where (a) is due to the fact that P(f2, p) is a concave function 
with respect to p. Thus, ^(O) is a convex set for any ft C V. 
It is easy to check that the sets S^x = {(P, P, p) : < p < 

Pmax}, ^2 - {(P,P,P) : T,Pi < Ph ^3 = {(R,Pp) ■ 

R < R}, and ^ 4 = {(P, P, p) : P < P tot } are also convex 
sets, and, as a result, /C = rincv\{<l}'^(^ U {s}) n ^1 n ^ 2 n 
^3 (~l £^4 is a convex set. ■ 
Having proved that © is a convex program, in order to 
use Proposition^ we need to check that the conditions of 
Definitions Q] |2] and |3] are satisfied. Part (i) of Definition 
Q] follows from the linearity of the objective in (0. For 
part (ii) of Definition [T] we specify an infeasibility measure 
Infeas(-) : Rl v l+ 2 -> M as follows^: 

Infeas((P, P, p)) = max jo, -p, p - p max , Rq-R,P- P to t 

E^-P,^- n min P(OU{ S },p)} 

(7) 

The conditions of part (ii) of Definition [T] are verified in the 
following theorem. 

Theorem 4. For a given vector (P, P, p) G M' v ' +2 and any 
e > we can either (a) determine in polynomial time if 
Infeas((R, P, p)) ^ £ and if not (b) find in polynomial time a 
vector c G Rl v l +2 , such that for every (P',P',p') satisfying 
Infeas((R',P',p')) e, c T (P',P',p') < c T (P,P,p). 

Proof: Part (a) requires checking that each of the ar- 
guments of the max of © is smaller than or equal to e in 
polynomial time. The first six terms are linear functions and 
can be easily computed. The last term can be compared to 
e by performing a minimization of a submodular function, 
which as was shown in Section IIV-BI can also be computed 
in polynomial time. 

We focus on condition (b). In this case Infeas((P, P, p)) > 
e, meaning that at least one of the arguments of the max of 
dT) is larger than e. We consider each case separately. 

If — pi > e we set c = — e^ +2 where has a one in the 
i th position and zeros everywhere else, which can be easily 
checked to satisfy the condition of part (b). 

Similarly, for the cases pi — p maXy i > £, Po — R > £, 
P - Ptot > e and Y^iPi - P we set c = e 4+2 , c = -ei, 
c = e 2 , c = (0, — 1, 1, . . . , 1) respectively. 

For the last case, let il* = argminfjgvUd} F(Sl U {s}, p). 
We have R-F(n*U{s}, p) > e. Since the function F(Q*, p) 

6 With a slight abuse of notation we use max on vector quantities by taking 
the maximum over the components of the vector. 
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is continuous and differentiable with respect to p, and the set 
<&(Sl*,e) = {(R,P,p) : R sC F(Q*,p)+e} is convex (which 
can be shown as in the proof of Lemma [TJ, then the vector 
c = (1, 0, — V p F(f2*, p)) is the normal to a hyperplane that 
separates (R,P,p) from the set ^(f2*,e). In other words, 
for all (R',P',p') e "(f(fi*) we have c T (R',P',p') < 
c T (R, P, p). Noting that {(R f , P' , p') : Infeas((i?', P', p')) < 
£} C ^(f2*,e), we conclude that part (b) holds in this case 
as well. ■ 
Having proved these preliminary results, we are ready to 
prove the main result of this section. 

Theorem 5. The optimization in (|6| can be solved in polyno- 
mial time on the size of the problem. 

Proof: The proof uses Proposition^ which requires 
verifying the convexity of the problem together with the 
conditions of polynomial computability, polynomial growth, 
and polynomial boundedness of the feasible set. Convexity 
was proved in Lemma [TJ while polynomial computability was 
shown in Theorem|4] Polynomial growth follows from the fact 
that F(Cl, p) is the log of a polynomial on p with degree at 
most |V|, while the objective and remaining terms that define 
the infeasibility measure are linear on (R, P, p). Finally, to 
check that feasible set if polynomially bounded, we note that 
the feasible set is a subset of the hypercube 

{(R,P,p):0^(R,P,p)^(R 

max • Ptot, Pmax)} 

where R max = min^ eV \{d} F(Cl U {s}, p max ). It follows that 
the feasible set is contained in the Euclidean ball centered 
at the origin with radius || (Rmax, Ptot, Pmax) ||2, which can 
be easily checked to grow polynomially on the size of the 
problem. ■ 
The general optimization problem $6fy can be specialized to 
a power minimization with a minimum rate constraint, and 
to a rate maximization with a total power constraint. Both 
problems can be solved in polynomial time, as stated in the 
following corollaries to TheoremTS] 

Corollary 2. The following power minimization problem can 
be solved in polynomial time. 

|V| 

minimize > pi 
p 

subject to Rq < F(Q. U {s}, p) for all Q, C V\{d} 

pmax ■ (8) 
Proof: The corollary follows from Theorem|5]by setting 

/il = 0, /i 2 = 1, Ptot = J2i=l Pmax,i- ■ 

Corollary 3. The following rate maximization can be solved 
in polynomial time. 

maximize R 

R,p 

subject to < R ^ F(Cl U {s}, p) for all il C V\{d} 

|V| 

2J Pi ^ Plot 
i=l 

< P < Pmax (9) 



Proof: The corollary follows from Theorem|5]by setting 
Hi = —1,^2 = and R = 0. ■ 

VI. Algorithms and simulations 

In this section we study different algorithms and and provide 
simulation results regarding submodular function minimization 
and power allocation problems. In the first subsection we 
look at minimum norm algorithm for submodular function 
minimization and we show that this algorithm can find the 
approximation to capacity of layered Gaussian relay networks 
with more that 300 nodes in a couple of minutes. In the 
second subsection we propose a heuristic algorithm to find 
the optimum power allocation for Gaussian relay networks. 

A. submodular function minimization 

One approach to solve the submodular minimization prob- 
lem due to Lovasz is based on extension of the set function 
/ : 2 V ->• R to a convex function g : [0,1] |v| -)■ K that 
agrees with / on the vertices of the hypercube [0, l]' v ', with 
a guarantee that min^cv f(A) is equal to min x g(x) for 
x e [0,1]' V L In this section we assume the normalization 
/(0) = 0. 

The Lovasz extension g of any set function / can be defined 
as follows. For a given x g [0,l]' v ' order the elements of 
V such that x{v±) ^ ^(^2) ^ ■•• x(v n ), where x(vi) is 
the t>ith element of the vector x. Set Ao = 1 — x(v\), \ = 
x(vi) - x(v i+ i), X n = x(v n ), and 

n 

#( X ) = ^2 X if({ v l, v 2,---,Vi})- 

Define lg = 6 R" and l{ Vl ,v 2 ....,vi} as an n dimensional 
vector such that the coordinates 1)2, ■ ■ ■ , v% are equal to one 
and all the other coordinates are equal to zero. Then, it is 

easy to see that x = Ya=o ^{ui,**,...,^} > E"=o A '' = 1 
and \ ^ 0. So, x is a unique linear convex combination 
of some vertices of the hypercube and g(x) is linear convex 
combination of values of / on those vertices. 

A key result is that if / is submodular its Lovasz extension 
g is a convex function |[T4l . ifTTl . In addition, finding the 
minimum of the submodular function / over subsets of V 
is equivalent to finding the minimum of the convex function 
g in the hypercube [0, 1]I V L The optimization can be done in 
polynomial time using Ellipsoid algorithm [14|. 

There are other algorithms with faster running time to solve 
the submodular minimization problem IfT^ll , IfTTl , |f20l . To 
the best of our knowledge, the running time of the fastest 
algorithm is in the order of 0(n 5 a + n 6 ), where a is the 
time that the algorithms takes to compute f(A) for any subset 
AC V [20|. For ADT networks, Gaussian relay networks, and 
erasure networks, a is the time to compute: the rank of n x n 
matrices, the determinant of n x n matrices, and equation (0), 
respectively. 

However, for networks of large size, a complexity of 
0(n 5 a + n 6 ) may still be computationally cumbersome. 
As a result, in these cases it is desirable to have faster 
algorithms. Recently, Fujishing IfTTl . ifTZl showed that the 
minimization of any submodular function can be cast as a 
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number of nodes 

Fig. 2. Running time of minimum norm algorithm for a layered Gaussian 
relay network. Each layer has four nodes. 



minimum norm optimization over the base polytope of /, 

Bf = P f n {x | £ ieV *W = /0>)}, where 



r> def I 



VA C V : ^2 x(i) < f(A) 1 

iGA J 



and the corresponding minimum norm optimization is 



minimize 



l x l|2, 



subject to x £ By. 
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Letting x* be the solution of this minimization, the set A* = 
{vi : x*(vi) < 0} is the solution to ruling f{A). Whether 
the above optimization problem can be solved in polynomial 
time is an open problem. However empirical studies [12| 
have shown that this algorithm has comparable or even faster 
running times than the other algorithms with polynomial time 
performance guarantees. 

In our specific setting, for layered Gaussian relay networks 
of size up to around 300 nodes with 4 nodes per layer, we 
were able to find the approximate capacity (cf. (O) in order 
of minutes (see Figure [2j. In order to solve the minimization 
( [Tol l we used the Matlab package provided in fl9l . 

B. Power allocation 

In Section [V] we have shown that the Ellipsoid method can 
be used to solve the optimization in (J8J in polynomial time. 
While in theory this result shows that the optimization in ([8J 
is tractable, in practice the Ellipsoid method has a number 
of shortcomings that limit its usability. On the one hand the 
running time of the Ellipsoid method can be large compared 
to alternative algorithms, and on the other hand, for high 
dimensional problems it has shown numerical instability. In 
this section, we propose a heuristic algorithm to solve the 
optimization in (J8j and show that this algorithm converges to 
the right solution. While the running time can be exponential 
on the network size, we show through simulations that the 
algorithm often converges within a time proportional to the 
network size. 



Our proposed algorithm, Algorithm Q] (see pseudo-code 
below), is based on the cutting plane methods [24] in convex 
optimization. The optimization in (JH]) contains exponentially 
many constraints of the form 

Ro < F(flU{s},p) for ft C V\{d}. (11) 

In Algorithm [T] we first find the min-cut corresponding to 
assigning maximum power to all nodes in the network: 

fii = ai-gmin ocv ^ d -jF(il U {s}, 

Then, we modify the optimization in (J8j by replacing the 
constraint (fTTT i with 

Ro < F(fii U{s},p). 

The resulting convex program can be easily solved since it 
contains few constraints. 

After optimization, we let p* be the optimum power allo- 
cation for the current set of constraints and set 



flj = argmin 



fiCV\{d} 



F(nu{s}, P *). 



We iteratively add the constraint 

Ro < F{QiU{s},p) 

to our set of constraints, and solve the optimization again. We 
stop if the new constraint is already in the set of constraints. 

Algorithm 1 Power minimization 

Input: Channel gain matrix H, desired rate R, vector of 
nodes' power constraints p max - 

Output: Min-cut, Power assignment p* that achieves approx- 
imate to rate R with minimum sum of powers. 

<£ <- {}, p* <- 

ft* <- min^cv\{d} F(Q U {s}, p max ) 
if R^F(n*U{s}, Pm ^) then 

while n* £ ^ and F(Q* U {s}, p*) < R do 

p* <- min Pi 
subject to: 
R sC F(n U {s}, p) for all ft € 

< P < Pmax 

n* <- min^cvud} F(Q U {s}, p*) 
end while 

return fl* U {s}, p* 
else 

print The constraints are infeasible. 
end if 

Since in each iteration the algorithm adds a new constraint 
to the constraint set and the number of constraints in (fTTI) is 
finite, the algorithm is guaranteed to find the optimum power 
in a finite number of iterations, which can be exponential. 

We used simulations to test the performance of the algo- 
rithm for networks of varying size n ranging from 10 to 40. 
For each n, we generated 300 random networks with channels 
gains drawn i.i.d. using a Af(0, 1) distribution. We set the 
desired transmission rate R = 4 in ([8j and set a maximum 
power constraint in each node to p max = 100. The results are 
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Fig. 3. Number of constraints in AlgorithmQ]when the algorithm terminates. Fi §- 5 - Minimum sum ° f P° wers for optimization GD when R = 4 and the 
The error bars represent one standard deviation. network ls S enerated randomly as described in the paper. 




number of nodes 

Fig. 4. Running time in seconds of the power optimization in Algorithm[T] 
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Fig. 6. Number of nodes in the simplified network. 



shown in Figures where the the vertical bars represent 
±1 standard deviation around the mean computed over the 
300 random networks. 

Figure [3] shows the number of iterations of Algorithm 1 as 
a function of the number of nodes in the network. We see that 
the number of iterations grows 0(n 3 / 2 ) with the number of 
nodes. 

In Figure [4] we present the simulation time as a function of 
the network size. We observe a running time that grows slower 
than 0(n e ). The figure also shows that the power optimization 
of networks of 40 nodes completes in less than two hours on 
an Intel Xeon quad core CPU running at 2.33 GHz. 

The minimum sum of powers to approximately achieve 
i?n = 4 for networks of different size is presented in Figure [5] 
with a blue solid line. Interestingly, the plot shows that the 
minimum power concentrates around the mean. In addition, 
the figure shows diminishing returns in total power savings 
resulting from increasing the network size. 



For the special case of the diamond relay network with TV 
relays, l26l shows that a fraction k/(k + 1) of the network 
capacity can be approximately achieved by using only k relays. 
In our setting of a general Gaussian relay network and sum- 
power minimization, we are interested in investigating whether 
it is possible to remove a large fraction of the nodes from the 
network without significantly affecting its performance. 

In order to determine which nodes to remove from the 
network, we solve the sum-power minimization and compare 
the optimal power allocation p* of each node i to a threshold 
Pth- All relays with p* < P t h are removed from the network. 
Let jV be the set of nodes with p* ^ P t h- We optimize the 
power allocation for the network with node set {sJU^/KUjd} 
and determine whether the problem is feasible for R = 4. If 
the problem is infeasible, we enlarge jV by adding more relay 
nodes in decreasing order of p*, until the problem becomes 
feasible. 

Figure [5] shows in the red dashed curve the resulting 
minimum sum of powers obtained by setting P t h = 1- Figure[6] 
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shows the corresponding size of the set jV . We observe in 
Figure [6] that the number of nodes with allocated power p* 
exceeding P t h = 1 (possibly including more relays to make 
the problem feasible) remains fairly constant as the size of 
the network n increases. This means that most of the power 
is allocated to a small subset of the nodes. 

Removing the remaining nodes from the network and opti- 
mizing the power allocation again over the resulting simplified 
network results in the minimum total power plotted in Fig- 
ure [3] This figure shows that even though the number of nodes 
in the simplified network remains approximately constant as 
n increases, the total power required to approximately achieve 
Ro = 4 decreases with n. This is due to the fact that larger 
n allows to choose the best relays for the simplified network. 
There is some performance loss in terms of total power due 
to network simplification but this loss may be compensated 
by power savings arising from turning off some of the relays. 
While we have not modeled the power consumption of the 
clocks, CPU and other subsystems required to keep a relay 
active, in practice they may become comparable to the power 
consumed by radio transmissions. This makes network sim- 
plification very useful in practice. 
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